function [V lambda] = PCA(allGPAShape)
estShape = getProcrusteMeanShape(allGPAShape);
[a n N] = size(allGPAShape);
D = zeros(n*2,N);
for i=1:N
    D(:,i) = [allGPAShape(1,:,i)-estShape(1,:) allGPAShape(2,:,i)-estShape(2,:)]';
end
S_trick = 1/(N-1)*(D'*D);
[V_trick, lambda] = eig(S_trick);
lambda = (diag(lambda));
V = zeros(2*n,N);
for i=1:N
    V(:,i) = D*V_trick(:,i)/norm(D*V_trick(:,i));
end
